Методи розв`язання крайових задач в тому числі жорстких краєвих завдань

Методи розв'язання крайових задач, у тому числі «жорстких» крайових задач

Методи Олексія Юрійовича Виноградова

1 Введення

На прикладі системи диференціальних рівнянь циліндричної оболонки ракети - системи звичайних диференціальних рівнянь 8-го порядку (після поділу приватних похідних).

Система лінійних звичайних диференціальних рівнянь має вигляд:

Y (X) = A (x) ∙ Y (x) + F (x),

де Y (x) - шукана вектор-функція задачі розмірності 8х1, Y (X) - похідна шуканої вектор-функції розмірності 8х1, A (x) - квадратна матриця коефіцієнтів диференціального рівняння розмірності 8х8, F (x) - вектор-функція зовнішнього впливу на систему розмірності 8х1.

Тут і далі вектора позначаємо жирним шрифтом замість рисочок над буквами

Крайові умови мають вигляд:

U ∙ Y (0) = u,

V ∙ Y (1) = v,

де

Y (0) - значення шуканої вектор-функції на лівому краї х = 0 розмірності 8х1, U - прямокутна горизонтальна матриця коефіцієнтів крайових умов лівого краю розмірності 4х8, u - вектор зовнішніх впливів на лівий край розмірності 4х1,

Y (1) - значення шуканої вектор-функції на правому краї х = 1 розмірності 8х1, V - прямокутна горизонтальна матриця коефіцієнтів крайових умов правого краю розмірності 4х8, v - вектор зовнішніх впливів на правий край розмірності 4х1.

У випадку, коли система диференціальних рівнянь має матрицю з постійними коефіцієнтами A = const, рішення задачі Коші має вигляд [Гантмахер]:

Y (x) = e ∙ Y (x ) + E ∙ e ∙ F (t) dt,

де

e = E + A (xx ) + A (Xx ) / 2! + A (X - x ) / 3! + ...,

де E це одинична матриця.

Матрична експонента ще може називатися матрицею Коші або матріціантом і може позначатися у вигляді:

K (x ← x ) = K (x - x ) = E .

Тоді рішення задачі Коші може бути записано у вигляді:

Y (x) = K (x ← x ) ∙ Y (x ) + Y * (x ← x ),

де Y * (x ← x ) = E ∙ e ∙ F (t) dt це вектор приватного рішення неоднорідної системи диференціальних рівнянь.

2 випадок змінних коефіцієнтів

Цей варіант розгляду змінних коефіцієнтів перевірявся в кандидатській дисертації.

З теорії матриць [Гантмахер] відомо властивість перемножаемості матричних експонент (матриць Коші):

e = E ∙ e ∙ ... ∙ e ∙ e ,

K (x ← x ) = K (x ← x ) ∙ K (x ← x ) ∙ ... ∙ K (x ← x ) ∙ K (x ← x ).

У випадку, коли система диференціальних рівнянь має матрицю із змінними коефіцієнтами A = A (x), рішення задачі Коші пропонується шукати за допомогою властивості перемножаемості матриць Коші. Тобто інтервал інтегрування розбивається на малі ділянки і на малих ділянках матриці Коші наближено обчислюються за формулою для постійної матриці в експоненті. А потім матриці Коші, обчислені на малих ділянках, перемножуються:

K (x ← x ) = K (x ← x ) ∙ K (x ← x ) ∙ ... ∙ K (x ← x ) ∙ K (x ← x ),

де матриці Коші наближено визначаються за формулою:

K (x ← x ) = E , Де Δ x = X - X .

3 Формула для обчислення вектора приватного рішення неоднорідної системи диференціальних рівнянь

Ця дуже проста формула ще не обраховано на комп'ютерах. Замість неї обраховувалися значно раніше виведена і набагато складніша формула, наведена в:

Чисельний метод перенесення крайових умов для жорстких диференціальних рівнянь будівельної механіки Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стор 1409-003r.pdf

Замість формули для обчислення вектора приватного рішення неоднорідної системи диференціальних рівнянь у вигляді [Гантмахер]:

Y * (x ← x ) = E ∙ e ∙ F (t) dt

пропонується використовувати наступну формулу для кожної окремої ділянки інтервалу інтегрування і тоді вектор приватного рішення на всьому інтервалі буде складатися з векторів, обчислених за формулою:

Y * (x ← x ) = Y * (x - X ) = K (x - X ) ∙ K (x - T) ∙ F (t) dt.

Правильність наведеної формули підтверджується наступним:

Y * (x - X ) = E ∙ e ∙ F (t) dt,

Y * (x - X ) = e ∙ e ∙ F (t) dt,

Y * (x - X ) = e ∙ F (t) dt,

Y * (x - X ) = E ∙ e ∙ F (t) dt,

Y * (x ← x ) = E ∙ e ∙ F (t) dt,

що потрібно було підтвердити.

Обчислення вектора приватного рішення системи диференціальних рівнянь проводитися за допомогою подання матриці Коші під знаком інтеграла у вигляді ряду та інтегрування цього ряду поелементно:

Y * (x ← x ) = Y * (x - X ) = K (x - X ) ∙ K (x - T) ∙ F (t) dt =

= K (x - X ) ∙ (E + A (x - T) + A (X - T) / 2! + ...) ∙ F (t) dt =

= K (x - X ) ∙ (E F (t) dt + A ∙ (X - T) ∙ F (t) dt + A / 2! ∙ (X - T) ∙ F (t) dt + ...).

Ця формула справедлива для випадку системи диференціальних рівнянь з постійною матрицею коефіцієнтів A = const.

Для випадку змінних коефіцієнтів A = A (x) можна використовувати прийом поділу ділянки (x - X ) Інтервалу інтегрування на малі подучасткі, де на подучастках коефіцієнти можна вважати постійними A (x ) = Const і тоді вектор приватного рішення неоднорідної системи диференціальних рівнянь Y * (x ← x ) Буде на ділянці складатися з відповідних векторів подучастков, на яких матриці Коші наближено обчислюються за допомогою формул з постійними матрицями в експонентів.

4 Метод «перенесення крайових умов» в довільну точку інтервалу інтегрування

Метод обраховано на комп'ютерах. По ньому вже зроблено 3 кандидатських фіз-мат дисертації.

Метод підходить для будь-яких крайових задач. А для "жорстких" крайових задач показано, що метод вважає швидше, ніж метод С. К. Годунова до 2-х порядків (в 100 разів), а для деяких «жорстких» крайових задач не вимагає ортонормірованія зовсім. Дивись:

Чисельний метод перенесення крайових умов для жорстких диференціальних рівнянь будівельної механіки
Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стор 1409-003r.pdf

Повне рішення системи диференціальних рівнянь має вигляд

Y (x) = K (x ← x ) ∙ Y (x ) + Y * (x ← x ).

Або можна записати:

Y (0) = K (0 ← x ) ∙ Y (x ) + Y * (0 ← x ).

Підставляємо це вираз для Y (0) в крайові умови лівого краю і отримуємо:

U ∙ Y (0) = u,

U ∙ [K (0 ← x ) ∙ Y (x ) + Y * (0 ← x )] = U,

[U ∙ K (0 ← x )] ∙ Y (x ) = U - U ∙ Y * (0 ← x ).

Або отримуємо крайові умови, перенесені в точку x :

U ∙ Y (x ) = U ,

де U = [U ∙ K (0 ← x )] І u = U - U ∙ Y * (0 ← x ).

Далі запишемо аналогічно

Y (x ) = K (x ← x ) ∙ Y (x ) + Y * (x ← x )

І підставимо цей вираз для Y (x ) В перенесені крайові умови точки x

U ∙ Y (x ) = U ,

U ∙ [K (x ← x ) ∙ Y (x ) + Y * (x ← x )] = U ,

[U ∙ K (x ← x )] ∙ Y (x ) = U - U ∙ Y * (x ← x ),

Або отримуємо крайові умови, перенесені в точку x :

U ∙ Y (x ) = U ,

де U = [U ∙ K (x ← x )] І u = U - U ∙ Y * (x ← x ).

І так в точку x переносимо матричне крайове умова з лівого краю і таким же чином переносимо матричне крайове умова з правого краю і отримуємо:

U ∙ Y (x ) = U ,

V ∙ Y (x ) = V .

З цих двох матричних рівнянь з прямокутними горизонтальними матрицями коефіцієнтів очевидно отримуємо одну систему лінійних алгебраїчних рівнянь з квадратною матрицею коефіцієнтів:

∙ Y (x ) = .

А у разі «твердих» диференціальних рівнянь пропонується застосовувати порядкове ортонормірованіе матричних крайових умов у процесі їх переносу в розглянуту точку. Для цього формули ортонормірованія систем лінійних алгебраїчних рівнянь можна взяти в [Березін, Жидков].

Тобто, отримавши

U ∙ Y (x ) = U ,

застосовуємо до цієї групи лінійних алгебраїчних рівнянь порядкове ортонормірованіе і отримуємо еквівалентну матричне крайове умова:

U ∙ Y (x ) = U .

І тепер вже в цей проортонормірованное порядково рівняння підставляємо

Y (x ) = K (x ← x ) ∙ Y (x ) + Y * (x ← x ).

І отримуємо

U ∙ [K (x ← x ) ∙ Y (x ) + Y * (x ← x )] = U ,

[U ∙ K (x ← x )] ∙ Y (x ) = U - U ∙ Y * (x ← x ),

Або отримуємо крайові умови, перенесені в точку x :

U ∙ Y (x ) = U ,

де U = [U ∙ K (x ← x )] І u = U - U ∙ Y * (x ← x ).

Тепер вже до цієї групи лінійних алгебраїчних рівнянь застосовуємо порядкове ортонормірованіе і отримуємо еквівалентну матричне крайове умова:

U ∙ Y (x ) = U .

І так далі.

І аналогічно поступаємо з посередником матричними крайовими умовами, які переносяться з правого краю у розглянуту точку.

В результаті отримуємо систему лінійних алгебраїчних рівнянь з квадратною матрицею коефіцієнтів, що складається з двох незалежно один від одного поетапно проортонормірованних матричних крайових умов, що вирішується будь-яким відомим методом для отримання рішення Y (x ) В аналізованої точці x :

∙ Y (x ) = .

5 Другий варіант методу «перенесення крайових умов» в довільну точку інтервалу інтегрування

Цей варіант методу ще не обраховано на комп'ютерах.

Запропоновано виконувати інтегрування за формулами теорії матриць [Гантмахер] відразу від деякої внутрішньої точки інтервалу інтегрування до країв:

Y (0) = K (0 ← x) ∙ Y (x) + Y * (0 ← x),

Y (1) = K (1 ← x) ∙ Y (x) + Y * (1 ← x).

Підставимо ці формули в крайові умови і отримаємо:

U ∙ Y (0) = u,

U ∙ [K (0 ← x) ∙ Y (x) + Y * (0 ← x)] = u,

[U ∙ K (0 ← x)] ∙ Y (x) = u - U ∙ Y * (0 ← x).

V ∙ Y (1) = v,

V ∙ [K (1 ← x) ∙ Y (x) + Y * (1 ← x)] = v,

[V ∙ K (1 ← x)] ∙ Y (x) = v - V ∙ Y * (1 ← x).

Тобто отримуємо два матричних рівняння крайових умов, перенесені у розглянуту точку x:

[U ∙ K (0 ← x)] ∙ Y (x) = u - U ∙ Y * (0 ← x),

[V ∙ K (1 ← x)] ∙ Y (x) = v - V ∙ Y * (1 ← x).

Ці рівняння аналогічно об'єднуються в одну систему лінійних алгебраїчних рівнянь з квадратною матрицею коефіцієнтів для знаходження рішення Y (x) у будь розглянутій точці x:

∙ Y (x) = .

У разі «твердих» диференціальних рівнянь пропонується наступний алгоритм.

Використовуємо властивість перемножаемості матриць Коші:

K (x ← x) = K (x ← x ) ∙ K (x ← x ) ∙ ... ∙ K (x ← x ) ∙ K (x ← x)

і запишемо вирази для матриць Коші, наприклад, у вигляді:

K (0 ← x) = K (0 ← x ) ∙ K (x ← x ) ∙ K (x ← x),

K (1 ← x) = K (1 ← x ) ∙ K (x ← x ) ∙ K (x ← x ) ∙ K (x ← x),

Тоді перенесені крайові умови можна записати у вигляді:

[U ∙ K (0 ← x ) ∙ K (x ← x ) ∙ K (x ← x)] ∙ Y (x) = u - U ∙ Y * (0 ← x),

[V ∙ K (1 ← x ) ∙ K (x ← x ) ∙ K (x ← x ) ∙ K (x ← x)] ∙ Y (x) = v - V ∙ Y * (1 ← x)

або у вигляді:

[U ∙ K (0 ← x ) ∙ K (x ← x ) ∙ K (x ← x)] ∙ Y (x) = u *,

[V ∙ K (1 ← x ) ∙ K (x ← x ) ∙ K (x ← x ) ∙ K (x ← x)] ∙ Y (x) = v *.

Тоді розглянемо ліве перенесене крайове умова:

[U ∙ K (0 ← x ) ∙ K (x ← x ) ∙ K (x ← x)] ∙ Y (x) = u *,

[U ∙ K (0 ← x )] ∙ {K (x ← x ) ∙ K (x ← x) ∙ Y (x)} = u *,

[Матриця] ∙ {вектор} = вектор.

Цю групу лінійних алгебраїчних рівнянь можна піддати порядковому ортонормірованію, яке зробить строчки [матриці] ортонормированном, {вектор} торкнуться не буде, а вектор одержить перетворення. Тобто отримаємо:

[U ∙ K (0 ← x )] ∙ {K (x ← x ) ∙ K (x ← x) ∙ Y (x)} = u * .

Далі послідовно можна записати:

[[U ∙ K (0 ← x )] ∙ K (x ← x )] ∙ {K (x ← x) ∙ Y (x)} = u * ,

[Матриця] ∙ {вектор} = вектор.

Аналогічно і цю групу лінійних алгебраїчних рівнянь можна піддати порядковому ортонормірованію, яке зробить строчки [матриці] ортонормированном, {вектор} торкнуться не буде, а вектор одержить перетворення. Тобто отримаємо:

[[U ∙ K (0 ← x )] ∙ K (x ← x )] ∙ {K (x ← x) ∙ Y (x)} = u * ,

Далі аналогічно можна записати:

[[[U ∙ K (0 ← x )] ∙ K (x ← x )] ∙ K (x ← x)] ∙ {Y (x)} = u * ,

[Матриця] ∙ {вектор} = вектор.

[[[U ∙ K (0 ← x )] ∙ K (x ← x )] ∙ K (x ← x)] ∙ Y (x) = u * .

Аналогічно можна проортонорміровать матричне рівняння крайових умов і для правого краю незалежно від лівого краю.

Далі проортонормірованние рівняння крайових умов:

[U ∙ K (0 ← x)] ∙ Y (x) = u * ,

[V ∙ K (1 ← x)] ∙ Y (x) = v *

як і раніше об'єднуються в одну звичайну систему лінійних алгебраїчних рівнянь з квадратною матрицею коефіцієнтів для знаходження шуканого вектора Y (x):

∙ Y (x) = .

6 методів додаткових крайових умов

Цей метод ще не обраховано на комп'ютерах.

Запишемо на лівому краї ще одне рівняння крайових умов:

M ∙ Y (0) = m.

Як рядків матриці M можна взяти ті крайові умови, тобто вираження тих фізичних параметрів, які не входять в параметри крайових умов лівого краю L або лінійно незалежні з ними. Це цілком можливо, оскільки у крайових задач стільки незалежних фізичних параметрів яка розмірність завдання, а в параметри крайових умов входить тільки половина фізичних параметрів задачі. Тобто, наприклад, якщо розглядається задача про оболонці ракети, то на лівому краї можуть бути задані 4 переміщення. Тоді для матриці М можна взяти параметри сил і моментів, яких теж 4, так як повна розмірність такого завдання - 8. Вектор m правій частині невідомий і його треба знайти і тоді можна вважати, що крайова задача вирішена, тобто зведено до задачі Коші, тобто знайдений вектор Y (0) з виразу:

∙ Y (0) = ,

тобто вектор Y (0) знаходиться з рішення системи лінійних алгебраїчних рівнянь з квадратною невиродженої матрицею коефіцієнтів, що складається з блоків U і M.

Аналогічно запишемо на правому краї ще одне рівняння крайових умов:

N ∙ Y (0) = n,

де матриця N записується з тих же міркувань додаткових лінійно незалежних параметрів на правому краї, а вектор n невідомий.

Для правого краю теж справедлива відповідна система рівнянь:

∙ Y (1) = .

Запишемо Y (1) = K (1 ← 0) ∙ Y (0) + Y * (1 ← 0) і підставимо в останню систему лінійних алгебраїчних рівнянь:

∙ [K (1 ← 0) ∙ Y (0) + Y * (1 ← 0)] = ,

∙ K (1 ← 0) ∙ Y (0) = - ∙ Y * (1 ← 0),

∙ K (1 ← 0) ∙ Y (0) = ,

∙ K (1 ← 0) ∙ Y (0) = .

Запишемо вектор Y (0) через зворотну матрицю:

Y (0) = ∙

і підставимо в попередню формулу:

∙ K (1 ← 0) ∙ ∙ = .

Таким чином, ми отримали систему рівнянь виду:

В ∙ = ,

де матриця У відома, вектори u і s відомі, а вектори m і t невідомі.

Розіб'ємо матрицю В на природні для нашого випадку 4 блоки і отримаємо:

∙ = ,

звідки можемо записати, що

В11 ∙ u + B 12 ∙ m = s,

B 21 ∙ u + B 22 ∙ m = t.

Отже, шуканий вектор m обчислюється за формулою:

m = B 12 ∙ (s - B11 ∙ u).

А шуканий вектор n обчислюється через вектор t:

t = B 21 ∙ u + B 22 ∙ m,

n = t + N ∙ Y * (1 ← 0).

У разі «твердих» диференціальних рівнянь пропонується виконувати почергове порядкове ортонормірованіе.

Запишемо наведену вище формулу

∙ K (1 ← 0) ∙ ∙ =

у вигляді:

∙ K (1 ← x 2) ∙ K (x 2 ← x 1) ∙ K (x 1 ← 0) ∙ ∙ = .

Цю формулу можна записати у вигляді поділу лівій частині на твір матриці на вектор:

[ ∙ K (1 ← x 2)] ∙ {K (x 2 ← x 1) ∙ K (x 1 ← 0) ∙ ∙ } =

[Матриця] ∙ {вектор} = вектор

[ ∙ K (1 ← x 2)] {K (x 2 ← x 1) ∙ K (x 1 ← 0) ∙ ∙ } =

Тут слід сказати, що подвектор t піддавати перетворенню не потрібно, так як неможливо, так як його первісне значення не відомо. Але подвектор t нам виявляється і не потрібен для виконання завдання.

Далі запишемо:

[[ ∙ K (1 ← x 2)] ∙ K (x 2 ← x 1)] {K (x 1 ← 0) ∙ ∙ } =

[Матриця] {вектор} = вектор

[[ ∙ K (1 ← x 2)] K (x 2 ← x 1)] {K (x 1 ← 0) ∙ } = .

І так далі.

В результаті почергового ортонормірованія отримаємо:

В ∙ = ,

∙ = .

Отже, шуканий вектор m обчислюється за формулою:

m = B 12 ∙ (s - B11 ∙ u).

7 Формула для початку рахунки методом прогонки С. К. Годунова

Ця формула обраховано на комп'ютерах в кандидатській дисертації.

Розглянемо проблему методу прогонки С. К. Годунова.

редположім, що розглядається оболонка ракети. Це тонкостінна труба. Тоді система лінійних звичайних диференціальних рівнянь буде 8-го порядку, матриця A (x) коефіцієнтів матиме розмірність 8х8, шукана вектор-функція Y (x) матиме розмірність 8х1, а матриці крайових умов будуть прямокутними горизонтальними розмірності 4х8.

Тоді в методі прогонки С. К. Годунова для такого завдання рішення шукається в наступному вигляді:

Y (x) = Y (X) c + Y (X) c + Y (X) c + Y (X) c + Y * (x),

або можна записати в матричному вигляді:

Y (x) = Y (X) ∙ c + Y * (x),

де вектори Y (X), Y (X), Y (X), Y (X) - це лінійно незалежні вектора-рішення однорідної системи диференціальних рівнянь, а вектор Y * (x) - це вектор приватного рішення неоднорідної системи диференціальних рівнянь.

Тут Y (X) = | | Y (X), Y (X), Y (X), Y (X) | | це матриця розмірності 8х4, а c це відповідний вектор розмірності 4х1із шуканих констант c , C , C , C .

Але взагалі то рішення для такої крайової задачі з розмірністю 8 (поза рамками методу прогонки С. К. Годунова) може складатися не з 4 лінійно незалежних векторів Y (X), а повністю з усіх 8 лінійно незалежних векторів-рішень однорідної системи диференціальних рівнянь:

Y (x) = Y (X) c + Y (X) c + Y (X) c + Y (X) c +

+ Y (X) c + Y (X) c + Y (X) c + Y (X) c + Y * (x),

І якраз труднощі і проблема методу прогонки С. К. Годунова і полягає в тому, що рішення шукається тільки з половиною можливих векторів і констант і проблема в тому, що таке рішення з половиною констант має відповідати умовам на лівому краї (стартовому для прогонки) при всіх можливих значеннях констант, щоб потім знайти ці константи з умов на правому краї.

Тобто в методі прогонки С. К. Годунова є проблема знаходження таких початкових значень Y (0), Y (0), Y (0), Y (0), Y * (0) векторів Y (X), Y (X), Y (X), Y (X), Y * (x), щоб можна було почати прогонку з лівого краю x = 0, тобто щоб задовольнялися умови U ∙ Y (0) = u на лівому краї за будь-яких значеннях констант c , C , C , C .

Зазвичай ця трудність «долається» тим, що диференціальні рівняння записуються не через функціонали, а через фізичні параметри і розглядаються найпростіші умови на найпростіші фізичні параметри, щоб початкові значення Y (0), Y (0), Y (0), Y (0), Y * (0) можна було вгадати. Тобто завдання зі складними крайовими умовами так вирішувати не можна: наприклад, завдання з пружними умовами на краях.

Нижче пропонується формула для початку обчислень методом прогонки С. К. Годунова.

Виконаємо порядкове ортонормірованіе матричного рівняння крайових умов на лівому краї:

U ∙ Y (0) = u,

де матриця U прямокутна і горизонтальна розмірності 4х8.

В результаті отримаємо еквівалентну рівняння крайових умов на лівому краї, але вже з прямокутною горизонтальної матрицею U розмірності 4х8, у якої будуть 4 ортонормированном рядки:

U ∙ Y (0) = u ,

де в результаті ортонормірованія вектор u перетворений у вектор u .

Як виконувати порядкове ортонормірованіе систем лінійних алгебраїчних рівнянь можна подивитися в [Березін, Жидков].

Доповнимо прямокутну горизонтальну матрицю U до квадратної невиродженої матриці W:

W = ,

де матриця М розмірності 4х8 повинна добудовувати матрицю U до невиродженої квадратної матриці W розмірності 8х8.

Як рядків матриці М можна взяти ті крайові умови, тобто вираження тих фізичних параметрів, які не входять в параметри лівого краю або лінійно незалежні з ними. Це цілком можливо, оскільки у крайових задач стільки незалежних фізичних параметрів яка розмірність завдання, тобто в даному випадку їх 8 штук і якщо 4 задані на лівому краю, то ще 4 можна взяти з правого краю.

Завершимо ортонормірованіе побудованої матриці W, тобто виконаємо порядкове ортонормірованіе і отримаємо матрицю W розмірності 8х8 з ортонормированном рядками:

W = .

Можемо записати, що

Y (0) = (М ) Транспонована = М .

Тоді, підставивши у формулу методу прогонки С. К. Годунова, отримаємо:

Y (0) = Y (0) ∙ з + Y * (0)

або

Y (0) = М ∙ з + Y * (0).

Підставимо цю останню формулу в крайові умови U ∙ Y (0) = u і отримаємо:

U ∙ [М ∙ з + Y * (0)] = u .

Звідси отримуємо, що на лівому краї константи c вже не на що не впливають, тому що

U ∙ М = 0 і залишається тільки знайти Y * (0) з виразу:

U ∙ Y * (0) = u .

Але матриця U має розмірність 4х8 і її треба доповнити до квадратної невиродженої, щоб знайти вектор Y * (0) з рішення відповідної системи лінійних алгебраїчних рівнянь:

∙ Y * (0) = ,

де 0 - будь-який вектор, у тому числі вектор з нулів.

Звідси отримуємо за допомогою оберненої матриці:

Y * (0) = ∙ ,

Тоді підсумкова формула для початку обчислень методом прогонки С. К. Годунова має вигляд:

Y (0) = М ∙ з + ∙ .

8 Другий алгоритм для початку рахунки методом прогонки С. К. Годунова

Цей алгоритм обраховано на комп'ютерах в кандидатській дисертації.

Цей алгоритм вимагає доповнення матриці крайових умов U до квадратної невиродженої:

Початкові значення Y (0), Y (0), Y (0), Y (0), Y * (0) знаходяться з рішення наступних систем лінійних алгебраїчних рівнянь:

∙ Y * (0) = ,

∙ Y (0) = , Де i = , , , ,

де 0 - вектор з нулів розмірності 4х1.

9 Заміна методу чисельного інтегрування Рунге-Кутта в методі прогонки С. К. Годунова

Ця заміна формул Рунге-Кутта на формулу теорії матриць обраховано на комп'ютерах в кандидатській дисертації.

У методі С. К. Годунова як показано вище рішення шукається у вигляді:

Y (x) = Y (X) ∙ c + Y * (x).

На кожній конкретній ділянці методу прогонки С. К. Годунова між точками ортогоналізації можна замість методу Рунге-Кутта користуватися теорією матриць і виконувати розрахунок через матрицю Коші:

Y (X ) = K (x - X ) ∙ Y (X ).

Так виконувати обчислення швидше, особливо для диференціальних рівнянь з постійними коефіцієнтами.

І аналогічно через теорію матриць можна обчислювати і вектор Y * (x) приватного рішення неоднорідної системи диференціальних рівнянь. Або для цього вектора окремо можна використовувати метод Рунге-Кутта, тобто можна комбінувати теорію матриць і метод Рунге-Кутта.

10 Метод половини констант

Цей метод поки не обраховано на комп'ютерах.

Вище було показано, що рішення системи лінійних звичайних диференціальних рівнянь можна шукати у вигляді тільки з половиною можливих векторів і констант. Була наведена формула для початку обчислень:

Y (0) = М ∙ з + ∙ .

З теорії матриць відомо, що якщо матриця ортонормированном, то її зворотна матриця є її транспонована матриця. Тоді остання формула набуває вигляду:

Y (0) = М ∙ з + U ∙ u

або

Y (0) = U ∙ u + М ∙ з

або

Y (0) = ∙ ,

Таким чином записана в матричному вигляді формула для початку рахунки з лівого краю, коли на лівому краї задоволені крайові умови.

Далі запишемо V ∙ Y (1) = v і Y (1) = K (1 ← 0) ∙ Y (0) + Y * (1 ← 0) спільно:

V ∙ [K (1 ← 0) ∙ Y (0) + Y * (1 ← 0)] = v

V ∙ K (1 ← 0) ∙ Y (0) = v - V ∙ Y * (1 ← 0)

і підставимо в цю формулу вираз для Y (0):

V ∙ K (1 ← 0) ∙ ∙ = V - V ∙ Y * (1 ← 0).

V ∙ K (1 ← 0) ∙ ∙ = P.

Таким чином ми отримали вираз виду:

D ∙ = P,

де матриця D має розмірність 4х8 і може бути природно представлена у вигляді двох квадратних блоків розмірності 4х4:

∙ = P.

Тоді можемо записати:

D 1 ∙ u + D 2 ∙ c = p.

Звідси отримуємо, що:

c = D 2 ∙ (p - D 1 ∙ u )

Таким чином, шукані константи знайдені.

Далі показано як застосовувати цей метод для вирішення «жорстких» крайових задач.

Запишемо

V ∙ K (1 ← 0) ∙ ∙ = P.

спільно з K (1 ← 0) = K (1 ← x 2) ∙ K (x 2 ← x 1) ∙ K (x 1 ← 0) і отримаємо:

V ∙ K (1 ← x 2) ∙ K (x 2 ← x 1) ∙ K (x 1 ← 0) ∙ ∙ = P.

Цю систему лінійних алгебраїчних рівнянь можна представити у вигляді:

[V ∙ K (1 ← x 2)] ∙ {K (x 2 ← x 1) ∙ K (x 1 ← 0) ∙ ∙ } = P.

[Матриця] ∙ {вектор} = вектор

[V ∙ K (1 ← x 2)] {K (x 2 ← x 1) ∙ K (x 1 ← 0) ∙ } = P .

І так далі.

У результаті почергового вичленення матриць ліворуч з вектора і ортонормірованія отримаємо систему:

D ∙ = P ,

Звідси отримуємо, що:

c = D 2 ∙ (p - D 1 ∙ u )

Таким чином, шукані константи знайдені.

11 застосовують формули ортонормірованія

Ці формули обраховані в кандидатській дисертації.

Взято з: Березін І.С., Жидков Н.П. Методи обчислень, том II, Державне видавництво фізико-математичної літератури, Москва, 1962 р. 635 стор

Нехай дана система лінійних алгебраїчних рівнянь порядку n:

А = .

Тут над векторами поставимо рисочки замість їх позначення жирним шрифтом.

Будемо розглядати рядки матриці А системи як вектори:

= ( , , ..., ).

Ортонорміруем цю систему векторів.

Перше рівняння системи А = ділимо на .

При цьому отримаємо:

+ + ... + = , = ( , , ..., ),

де = , = , = 1.

Друге рівняння системи замінюється на:

+ + ... + = , = ( , , ..., ),

де = , = ,

= - ( , ) , = - ( , ) .

Аналогічно робимо далі. Рівняння з номером i прийме вигляд:

+ + ... + = , = ( , , ..., ),

де = , = ,

= - ( , ) - ( , ) - ... - ( , ) ,

= - ( , ) - ( , ) - ... - ( , ) .

Процес буде здійснимо, якщо система лінійних алгебраїчних рівнянь лінійно незалежна.

В результаті ми прийдемо до нової системи С = , Де матриця С буде з ортонормированном рядками, тобто має властивість С * З = E, де Е - це одинична матриця.

(Таким чином, рішення системи можна записати у вигляді = С .)

12 Висновок формул, запозичений з «Теорії матриць» Гантмахера

Система лінійних звичайних диференціальних рівнянь з постійними коефіцієнтами має вигляд:

Y (X) = A Y (x) + F (x). (1)

Розкладемо Y (x) в ряд Маклорена за ступенями x:

Y (x) = Y + Y x + Y x / 2! + ..., Де Y = Y (0), Y = Y (0), ... (2).

З (1) почленно диференціюванням при А = const і F (x) = 0 отримаємо:

Y = AY = A Y, Y = AY = A Y, (3)

Поклавши в (3) x = 0 і підставивши в (2) отримаємо:

Y (x) = Y + Ax Y + A x / 2! Y + ... = E Y , (4)

де e = E + Ax + A x / 2! + ..., Де Е - одинична матриця. (5)

Якщо взяти x = x , То (4) заміниться на

Y (x) = e Y (x ), (6)

Розглянемо випадок A = const і F ≠ 0.

Введемо в розгляд вектор-функцію Ya (x) у вигляді: Y (x) = e Ya (x). (7)

Продіффренціруем (7) і підставимо в (1). Отримаємо:

e Ya (X) = F (x). (8)

При отриманні (8) враховувалося, що:

= = A + A x + A x / 2! + ... = A e .

З (8) випливає, що:

Ya (x) = c + . (9)

Підставимо в (7) і отримуємо:

Y (x) = e c + e . (10)

Поклавши x = x в (10) отримаємо:

c = e Y (x ). (11)

Остаточно отримуємо:

Y (x) = e Y (x ) + E . (12)

Мій батько запропонував використовувати і іншу (набагато ефективнішу за часом рахунку) матричну формулу замість матричної експоненти - щось на основі Вольтерра. Це є в статті в журналі «Математичне моделювання»:

Я тут якось днями поговорив з батьком і ось що цікаве щодо методу з'ясувалося.

Коли я сам вирішував своїм методом «перенесення крайових умов» «жорсткі» крайові задачі, то я в кожній розглянутій точці x = x * вирішував відповідну систему лінійних алгебраїчних рівнянь для знаходження рішення Y (x *) в кожній розглянутій точці x = x *.

А мій батько стверджує, що «жорсткими» бувають тільки крайові задачі, а початкові завдання «жорсткими» не бувають. І тому він знаходив рішення Y (x *) в якійсь одній точці x = x * по моєму методу, а далі вирішував (вліво і вправо від даної точки x = x *) як завдання Коші (від знайдених початкових умов Y (x *) в цій одній точці): Y (x) = K (х ← x *) · Y (x *) + Y * (x ← x *). І якщо це так, то так вирішувати, природно, набагато швидше, оскільки треба тільки домножать на матрицю Коші (матріціант, матричну експоненту) замість вирішення систем лінійних алгебраїчних рівнянь.

Я думаю, що це дуже істотне уточнення по швидкості рахунку.

P. P. P. S. Метод рішення "жорстких" початкових завдань, тобто «жорстких» завдань Коші. Придуманий ввечері 16 січня 2008 року.

Метод перенесення крайових умов працює успішно. Отже, за аналогією можна намагатися вирішувати «жорсткі» початкові завдання, тобто «жорсткі» задачі Коші.

Нехай дана початкова задача:

Y (X) = A · Y (x), Y (0) = Y поч,

Можемо записати:

Y (0) = K (0 ← x) · Y (x),

K (0 ← x) · Y (x) = Yнач

K (0 ← x3) · K (x3 ← x2) · K (x2 ← x) · Y (x) = Yнач

[K (0 ← x3)] · {K (x3 ← x2) · K (x2 ← x) · Y (x)} = Yнач

[Матриця] {} вектор вектор

Виконуємо порядкове ортонормірованіе цієї системи лінійних алгебраїчних рівнянь з квадратною матрицею коефіцієнтів і отримуємо систему з ортонормированном рядками в квадратної матриці:

[K (0 ← x3)] орт · {K (x3 ← x2) · K (x2 ← x) · Y (x)} = Yнач_орт

Аналогічно записуємо

[[K (0 ← x 3)] орт · K (x 3 ← x 2)] · {K (x 2 ← x) · Y (x)} = Y поч

[Матриця] {} вектор вектор

Далі виконуємо порядкове ортонормірованіе і отримуємо:

[[K (0 ← x 3)] орт · K (x 3 ← x 2)] орт · {K (x 2 ← x) · Y (x)} = Y нач_2орт

Аналогічно отримуємо

[[K (0 ← x 3)] орт · K (x 3 ← x 2)] орт · K (x 2 ← x)] орт · Y (x) = Y нач_3орт

Тобто отримали підсумкову систему лінійних алгебраїчних рівнянь для знаходження вектора Y (x). Для цього вектор Yнач_3орт треба домножити зліва на матрицю транспоновану по відношенню до лівої ортонормованій матриці, тому що якщо матриця ортонормірованна, то зворотна їй є транспонована матриця.

P. P. P. P. S. 11 вересня 2009:

Довго було лінь записувати, але раптом хтось сам не відразу здогадається - ще один метод вирішення "жорстких" початкових завдань, тобто «жорстких» завдань Коші. Хоча у мене є така підозра, що «жорсткими» бувають тільки крайові задачі, а початкові завдання «жорсткими» НЕ бувають. Але про всяк випадок наведу метод рішення "жорстких" початкових завдань.

Початкові умови мають вигляд:

Y (0) = Y поч.

Повне рішення системи звичайних диференціальних рівнянь Y (X) = A · Y (x) + F (x) має вигляд:

Y (x) = K (x ← 0) · Y (0) + Y * (x ← 0).

Або можна записати:

Y (0) = K (0 ← x 1) · Y (x1) + Y * (0 ← x 1).

Підставимо це вираження в крайові умови і отримаємо:

K (0 ← x 1) · Y (x1) + Y * (0 ← x 1) = Y поч

або

K (0 ← x 1) · Y (x1) = Y поч - Y * (0 ← x 1)

або

K 1 · Y (x 1) = Y 1.

Проортонорміруем цей вислів порядково і отримаємо еквівалентну вираз:

K 1орто · Y (x 1) = Y 1орто.

Тоді

Y (x 1) = (K 1орто) транспонирует · Y 1орто.

Підставимо замість Y (x1) вираз через Y (x2) і отримаємо:

K (x 1 ← x 2) · Y (x2) + Y * (x 1 ← x 2) = (K 1орто) транспонирует · Y 1орто

або

K (x 1 ← x 2) · Y (x2) = (K 1орто) транспонирует · Y 1орто - Y * (x 1 ← x 2)

або

K 2 · Y (x 2) = Y 2.

Проортонорміруем порядково і отримаємо еквівалентну вираз:

K 2орто · Y (x 2) = Y 2орто.

Тоді:

Y (x 2) = (K 2орто) транспонирует · Y 2орто.

І так далі.

P. P. P. P. P. S. Метод для чисельного інтегрування диференціальних рівнянь.

Читали нам якось в Бауманки чисельні методи рішення диференціальних рівнянь. І, здається, приводили аналітичний висновок формул одного з авторів. Чи це просто майнуло в підручнику (я маю на увазі висновок формул). Вже не дуже пам'ятаю. Запам'яталася тільки власна думка, що людям взагалі-то найпростіше даються геометричні аналогії і висновки, зроблені на основі зрозумілих геометричних картинок. Ну, ось тоді я і намалював один з варіантів чисельного рішення диференціальних рівнянь і пам'ятаю навіть переклав геометричні малюнки в літерні формули наближених обчислень. Зараз повторно виводити літерні формули для чисельного інтегрування диференціальних рівнянь мені не здається цікавим. А ось привести картинки тих студентських думок цілком можна для обговорення.

Далі йде картинка з текстом і з малюнками чисельного інтегрування: